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TIME -ASYMPTOTIC SOLUTIONS OF THE 
NAVIER- STOKES EQUATIONS FOR FREE SHEAR FLOWS USING 
AN ALTERNATING- DIRECTION IMPLICIT METHOD 

David H. Rudy and Dana J. Morris 
Langley Research Center 

SUMMARY 

An uncoupled time-asymptotic alternating-direction implicit method for solving the 
Navier-Stokes equations was tested on two laminar parallel mixing flows. A constant 
total temperature was assumed in order to eliminate the need to solve the full energy 
equation; consequently, static temperature was evaluated by using an algebraic relation- 
ship. For the mixing of two supersonic streams at a Reynolds number of 10^, converged 
solutions were obtained for a time step 5 times the maximum allowable size for an explicit 
method. The solution diverged for a time step 10 times the explicit limit. Improved con- 
vergence was obtained when upwind differencing was used for convective terms. Larger 
time steps were not possible with either upwind differencing or the diagonally dominant 
scheme of Khosla and Rubin. Artificial viscosity was added to the continuity equation in 
order to eliminate divergence for the mixing of a subsonic stream with a supersonic stream 
at a Reynolds number of 10^. 


INTRODUCTION 

This report presents the results of an investigation using an uncoupled alternating- 
direction implicit (ADI) finite-difference technique to obtain steady-state solutions for two- 
dimensional, high Reynolds number, compressible free shear flows by using a time- 
asymptotic approach. Preliminary results from this study are presented in reference 1. 

The review papers of Peyret and Viviand (ref. 2) and Taylor (ref. 3) summarize the 
published finite-difference solutions of the viscous compressible time -dependent Navier- 
Stokes equations. Most of the published solutions have used explicit differencing schemes. 
Although they are conceptually less complex and thus more easily coded than implicit 
methods, explicit finite -difference methods are restricted to small time steps relative to 
the spatial grid size for numerical stability. Consequently, explicit methods generally 
require a large total computation time to reach a steady -state flow condition. On the other 
hand, most implicit methods are found to be unconditionally stable in linearized stability 
analyses for model equations similar to the Navier-Stokes equations. Therefore, it is 



expected that implicit methods will allow larger time Steps than explicit methods when 
applied to the Navier -Stokes equations. 

The ADI technique developed by Peaceman and Rachford (ref. 4) and Douglas (ref. 5) 
was used in the present study. This ADI method is a two-step procedure requiring reduc- 
tion of tridiagonal matrices for which an efficient solution algorithm, the Thomas algorithm 
(ref. 6), exists. The method was originally applied to the two-dimensional heat conduction 
equation in reference 4 and was later applied to a system of hyperbolic equations by 
Gourlay and Mitchell (ref. 7). For both of these model problems, the ADI method was 
shown to possess unconditional stability. Roache (ref. 8) states that ADI methods are cur- 
rently the most popular approach to computing incompressible viscous Navier-Stokes flows. 
He notes that although these methods have not been found to be unconditionally stable in 
practice, the somewhat larger time steps obtainable have resulted in faster overall com- 
putation times (by a factor of 2 or more) than explicit methods for such flows. 

The AlDI method has been applied to the time -dependent compressible Navier-Stokes 
equations in only a few papers. In 1966, Polezhaev (ref. 9), for example, used an ADI 
method to obtain solutions for a two-dimensional natural convection problem. His uncou- 
pled ADI procedure removed the diffusion time-step limitation; however, he found that a 
sufficient condition for stability limited the time step to the usual maximum explicit value. 
In 1973, Briley and McDonald (ref. 10) presented a coupled ADI method in which nonlin- 
earities at the implicit time level were linearized by using a Taylor's series expansion 
about the known time level. The resulting system of linear difference equations was 
solved simultaneously by using a Douglas-Gunn (ref. 11) ADI approach. The method was 
found to be stable for very large Courant numbers = 1 corresponds to the Courant- 

Freidrichs-Lewy condition j in the calculation of three-dimensional subsonic flow in a 
straight duct with a rectangular cross section. For a flow with a Mach number of 0.044 and 
a Reynolds number of 60, stable solutions were obtained for Courant numbers up to 1250. 
For a Mach number of 0.5 and a Reynolds number of 600, the time step was gradually 
increased as the solution progressed; this gave an average Courant number of 73. Thus, 
the maximum Courant number decreased with increasing Reynolds number, perhaps 
because of diagonal dominance problems in the coefficient matrix as suggested in refer- 
ence 10. The computational effort per time step was reported to be twice that of most 
explicit methods. Baum and Ndefo (ref. 12) have also published an implicit method based 
on the Peaceman -Rachford procedure. This method iteratively solves coupled nonlinear 
difference equations by using a quasi -linearization technique. Reference 12 does not 
include results from solutions of the compressible Navier-Stokes equations. However, in 
reference 13 the method has been applied to the laminar hypersonic near wake of a sharp 
cone. 
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The ADI method investigated in the present paper is an uncoupled procedure which 
is applied to the mixing of two parallel streams. The results of reference 1 for the mix- 
ing of two supersonic streams have been improved for central differencing and extended 
to include computations with two other forms of differencing which give unconditional 
diagonal dominance in the coefficient matrices. In addition, results are presented for the 
mixing of a subsonic stream with a supersonic stream. 

SYMBOLS 


H 

L 

M 

N, 


Co 


N 


Re 


P 

R 

S 

T 

t 


'^ref 


speed of sound 


enthalpy 


reference length (L = 2.54 cm) 


Mach number 
Courant number, 

Reynolds number. 


At 


Ax 


lul + Iv I + \/2c 
Ps^refL 




pressure 
gas constant 
= 198.6/Tg 
temperature 
time 

streamwise velocity (x-direction) 
reference velocity, ^2Hg 

free-stream velocity on low-velocity side of mixing region 
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U2 free-stream velocity on high-velocity side of mixing region 

V normal velocity (y -direction) 

x,y coordinates in streamwise and normal directions, respectively 

a artificial diffusion coefficient 

y ratio of specific heats 

p molecular viscosity 

p density 

maximum values of x and y, respectively, in solution domain 

rn«.x ni«-x 
Subscripts: 

i,j index denoting grid -point spatial location 

s stagnation condition 

t,x,y derivative with respect to time, x-direction, and y-direction, respectively 

Superscripts: 

n,* index denoting time level 

Abbreviations : 

ADI alternating-direction implicit 

IC initial condition 

A bar over a symbol denotes a dimensional quantity. 
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GOVERNING EQUATIONS 


The nonconservative form of the governing partial differential equations was used in 
the present investigation. These equations are written as follows: 

Continuity 


Pj. + pVy + pUjj + vpy + VLp^ = 0 


x-momentum 


pu. . pvuy . puu, = -Px ^ - 4; K)^ * 5 -^ [p(>V - v,)] 


( 1 ) 


( 2 ) 


y -momentum 


pv, . pwy . puv, = -p^ . " "4 


( 3 ) 


These equations are nondimensionalized with respect to a reference length and stagnation 
flow conditions, that is, 


= _P_ 


T = ^ 


P = 




x = 2 L 


'^ref ~ 


u = 


V = 


u 

u 


^ref 

V 

V 

l/2Hs 

uref 




N. 


Re 


Ps^ref^ 


Ps 
lu 


t = 


ref 


The pressure was evaluated by means of the perfect gas equation of state 
p = pRT 


(4) 

5 


II 


where R = 2Lz_i, Air (y = 1.4) was the test fluid. Only laminar (molecular) viscous 
2y 

effects were considered, with the Sutherland formula being used to express the viscosity 
as a function of temperature as follows : 


^ ^ ^3/2 liS 
T + S 


(5) 


Calculations for a Mach 3 jet into still air with the quasi-parallel code of Oh (ref. 14), 
which includes the energy equation, showed that there was less than a 5-percent variation 
in the total enthalpy throughout the mixing region from the constant value assumed in other 
calculations. For the conditions of case 1 (supersonic -supersonic parallel mixing), the 
calculations of Oh showed less than a 1-percent variation in total enthalpy. This small 
variation had a negligible effect on the other flow parameters. Thus, to simplify the sys- 
tem of governing equations and to reduce required machine storage, a constant total tem- 
perature of 294 K (530° R) was assumed. Therefore, the static temperature was evaluated 
by use of the algebraic relationship 

T = 1 - (u^ + v^) (6) 

which eliminated the need for solving the complete energy equation. Constant static pres- 
sure was assumed in order to generate initial values of density by using equations (4) 
and (6) along with the given initial velocities. The linearized version of equations (1) to (6) 
with the viscous terms neglected has been shown by Gottlieb and Gustafsson (ref. 15) to be 
well-posed for the initial value problem. Polezhaev (ref. 9) considered the nonconserva- 
tive form of the momentum and energy equations but used the conservative form of the con- 
tinuity equation. Briley and McDonald (ref. 10) used the conservative form except for the 
energy equation in their coupled procedure. 


DESCRIPTION OF NUMERICAL TECHNIQUE 


Solution Procedure 

In the ADI procedure used in the present investigation, a sequential (uncoupled) solu- 
tion of the finite -difference form of equations (1) to (3) is obtained for each interior row of 
grid points during the first one-haif time step (horizontal sweep) and for each interior col- 
umn of grid points during the second one-haif time step (vertical sweep). All spatial 
derivatives were approximated by centered (second-order) finite differences, and time 
derivatives were approximated by backward (first-order) differences. The nonlinear coef- 
ficients in the convective terms were lagged one-half time step. In addition, the pressure 
terms and cross -derivative terms were also treated explicitly in each sweep. At each 


6 



row or column the order of solution was as follows: (1) solve the x-momentum equation 
for u, (2) solve the y -momentum equation for v, and (3) solve the continuity equation for 
p. The temperature and viscosity were then updated for the entire flow field after each 
sweep. No iteration was performed at any time step during the calculation. After the 
horizontal sweep, the values of u, v, and p were updated for the two boundary rows of 
points (i.e., the first and last rows) according to the appropriate boundary conditions. 
Similarly, the first and last columns were updated after the vertical sweep. The solution 
was then marched in time, with steady state considered to be the asymptotic limit. This 
uncoupled technique was chosen to avoid the coding complexity and possible additional 
computing time associated with the block-tridiagonal structure of a fully coupled approach. 


Finite -Difference Forms of Equations 

x-momentum equation. - To illustrate the present ADI procedure, the finite -difference 
form of the x-momentum equation is shown. For the horizontal sweep, from time level n 
to an intermediate time denoted by *, with a uniform spatial grid, 





u, 


1+1,1 ~^i-l,.i ' ^ _ 

2 Ax / 


- p!" 
1 + 1 ,3 Iz 


Li 


2 Ax 


%e 

1 

NRe 

4 

3NRe 

2 


/p" . , + iJ} /u? . , - u" A /u? . + iJ}. / 

^1+1 b3 b]+l IJ _ IkJ LLili h2zl 

\ 2 /\ Ay / \ 2 /\ Ay y 


, n 


- v^ 




-y^ 


n [ '^i+1,1+1 ~ ''i-1,1 + 1 j n ( ''i+1,1-1 ~ ''i-1,1-1 

iJ+lV 2 Ax / i,j-l\ 2 Ax J 


- <i ) ( i^li * - vi.i ' 


3NRe Ax 


Ax 


n ^i+l,j+l ~ ^i+l,j-l^ _ n ( ^i-l,i+l ~ '^i-l,j-l 
^i+l,iV 2 Av / ^i-1,3 


2 / \ Ax 

-n 

2 Ay 


For the jth row, this equation has the tridiagonal matrix form 

A.u. H • + B.u. . + C.u. - . = D. 

1 1+1,3 1 


( 7 ) 


( 8 ) 
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where A^, B^, and Cj are matrix coefficients which are defined as follows: 




1 2 Ax 3N (Ax) 2 

Re 


lU 


(9a) 


4 

B. = + » 

^ At 3Nj^g(Ax)4 


f jj,? . + . A 

liJ i-l>3 + 


'* ^1+1,1 ^1,1 


3Npe(Ax)2 \ 2 


(9b) 




1 2 Ax 3Nj^g(Ax)2 


/m" . - + m" 
i+l>J ^ 


(9c) 


and D. is a vector of quantities which are evaluated at time n. The unknowns are 

Ui j Uj j, and j- Similarly, for the vertical sweep, from * to n + 1, the 

x-momentum equation in finite-difference form is given by 


n-H 1 * 

/u, T - U. 4’ 




_U iil 


b]\ At/2 

1 


* * 


1,1+1 i,]-i 

V 2 Ay y 


(' * \ / 5fc * \ 

^+l,j-Vl,A .fa+ M-Pi-i,j 
2 Ax / V 2 Ax J 


NRe Ay 
1 

NpeAy 




(10) 


For the ith column, this equation has the tridiagonal matrix form 


Vu-i * 


(H) 


where A^, B and Cj are matrix coefficients which are defined as follows : 


* * 

-p. .V. 

A. = b3 1,1 _ 


1 2 Ay Nj^g(Ay)2 


/ * * \ 


(12a) 
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« * 

2pi . 

B. = — + 
3 At 


NR^(Ay)' 


/♦ ♦\/=K * \ 

\i * \l-l ) 


* * 

c. = 1 


J 2 Ay Nj^g(Ay)^ 


' U. . + LL, . ^ ' 


(12b) 

(12c) 


and Dj is a vector of quantities which are evaluated at time *. The unknowns are 

u?t^, and The accuracy of all the difference equations is 

0(At,(Ax)2,(Ay)2). A more accurate transient solution was considered unnecessary since 
only the steady state was of interest. 

y-momentum equation. - For the jth row during the horizontal sweep, the 
y-momentum equation is written 


5k jjC 

AiVi.l,, + t = D. 


where 




A - -'’l.fl.i 1 


‘ 2 Ax Nn_(Ax)2 


! p? . + jU? . J 


'Re' 
1 


2p(‘i 

B. = — iii + - 

At Nr^(Ax)2 




.n 


D . .U . - 

c. = ^^3 ^>3 


1 K-1,3-^,3 


(13) 


(14a) 


(14b) 


(14c) 


^ 2 Ax Nj^g(Ax)2 \ 2 / 

For the ith column during the vertical sweep, the y-momentum equation has the form 


^fi,j-\ ^fi,i+i - 


where 


5k 5k 

= 

3 9. A\T -NT ^A,r\2 


u. . + a. . ^ 1 
hi + 


2 Ay Njjg(Ay)^ 


3NRe(^y)' 


/ * * \ 
Ki ^ ^3-1 


(15) 


(16a) 
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« * 

2pi . 1 

B. = — ^ 




3NRe(Ay)' 


* * \ / ♦ ♦ N 

\j ^ij-l ) ^ Pi,j ^i,j+l 


(16b) 


* * 

C ^ _ 1 


/ * * 
fM. i + P. 




i,j %j+l ) ^ 2 rij ^ ^ij+l , 

2 f 3Npg(Ay)2 \ 2 ) 


(16c) 


Continuity equation. - Since the continuity equation (eq. (1)) contains no viscous terms, 
it was found to be beneficial in some cases to add artificial viscosity terms to equation (1) 
for numerical stability. These artificial viscosity terms were incorporated explicitly so 
that the continuity equation becomes 


Pt + PUjj. + uPj^ + pVy + vPy = 0'|(^x)^p^ + (Ay)^Pyy 


(17) 


However, as described later, a = 0 was used in many calculations. The finite-difference 
form of equation (17) for the horizontal sweep is then given by 


Vi-l.) * = D; 


(18) 


where 


__Ll + 0 / 

2 Ax 

(19a) 

^ - 2q; 
At 

(19b) 

* 

Ai + a 

2 Ax 

(19c) 


For the vertical sweep, the continuity equation is given by 




(20) 
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where 


n+1 




2 Ay 


+ a 





n+1 

C. = IiJ- + a 

] 2 Ay 


(21a) 

(21b) 

(21c) 


TEST PROBLEMS 

Figure 1 shows the mixing problems (cases 1 and 2) chosen for use as the standard 
test problems. In these test problems, the flow is the mixing of two parallel streams 
where the computational region begins far downstream of the base of the infinitely thin 
splitter plates and is truncated farther downstream. Such a calculation obviously does 
not require the full Navier-Stokes equations, since solutions can be obtained with the con- 
ventional quasi-parallel approach, that is, using the boundary -layer equations with free 
shear flow boundary conditions. However, the test problems were chosen because solutions 
could be obtained with this alternate method for comparison with computed Navier-Stokes 
results. The test problems were also chosen because the parabolic solutions could be 
used to formulate and test boundary conditions. Calculations were made for the mixing of 
two supersonic (Mach 3 and 1.68) streams (case 1) as well as the mixing of a subsonic 
(Mach 0.11) stream and a supersonic (Mach 3) stream (case 2). Calculations for case 1 
were made for values of of 10^, 5.0 x 10^, and 8.1 x 10^. In all calculations for 

case 2, was 10^. 

BOUNDARY CONDITIONS 

A sketch of the computational domain indicating the boundary conditions which were 
specified for case 2 is presented in figure 2. For case 1, the appropriate supersonic 
boundary conditions were applied along the entire inflow and downstream outflow 
boundaries. 

The subsonic boundary conditions were chosen on the basis of the analysis given by 
Gottlieb and Gustafsson in reference 15. At the subsonic outflow boundary in case 2, the 
one-dimensional analysis indicated that one function value, either p or u, must be speci- 
fied. Results of tests with three different outflow boundary conditions obtained by using an 
explicit method (hopscotch) given in reference 1 showed that specifying the steady -state 
density was the best choice. Such a boundary condition is obviously not convenient for 
most applications since these downstream steady -state values are generally not known 
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a priori. In the present investigation, however, these values were obtained from the para- 
bolic code calculations. 

The inflow profiles of u and v were also taken from parabolic code calculations. 
The temperature was obtained from equation (6). The density was then computed by using 
equation (4) with the constant value of static pressure, p = 0.00389. 

Figure 2 also shows the arrangement of grid points at the upper and lower bound- 
aries. At each of these locations the boundary is positioned between two rows of grid 
points. The specified values of u and v at the upper inflow are again obtained from 
the parabolic code calculations. 

The linear extrapolation condition used at the downstream boundary is given for uni- 
form grid spacing by 

%I,j = %I-l,j ■ %I-2,j ^22) 

where f is either p, u, or v and NI is the value of i at the downstream boundary. 
This is equivalent to saying that the second derivative of f is zero at NI - 1. 

RESULTS AND DISCUSSION 

Supersonic-Supersonic Parallel Mixing (Case 1) 

Calculations were made for case 1 with Nj^g = 10^, 5.0 x 10^, and 8.1 x 10^. The 
four initial flow-field conditions that were used are as follows : 

IC 1 Parabolic code inflow values for u and v specified at inflow station and outflow 
station; interior values obtained by linear interpolation along each row; initial den- 
sity computed from equation (4) with assumption of constant static pressure; initial 
temperature computed from equation (6). 

IC 2 Parabolic code inflow values for u and v specified at inflow station; outflow- 
station parabolic code values of u and v multiplied by 0.8; interior values 
obtained by linear interpolation along each row; initial density computed from equa- 
tion (4) with assumption of constant static pressure; initial temperature computed 
from equation (6). 

IC 3 Same as IC 2 except outflow- station parabolic code values of u and v multiplied 
by 0.6. 

IC 4 Parabolic code inflow values for u and v specified at inflow station; initial den- 
sity computed from equation (4) with assumption of constant static pressure; all 
interior columns and the outflow column set equal to first column; initial tempera- 
ture computed from equation (6). 
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At 



1 



i 0.01 



(25) 


This convergence criterion is identical to that of reference 1. 

Calculations were presented in reference 1 for = 10^. As shown in figure 6 

of reference 1, the pressure profile exhibited pointwise oscillations. The inflow (x = 0) 
profile for that calculation was taken as the first computed profile in the parabolic solu- 
tion. At X = -0.075 in the parabolic solution, the initial u profile was selected to be 
an error function between the desired outer stream velocities, u^ and Ug. The initial 
V was assumed to be zero. Marching steps of Ax = 0.075 were taken. At the first com- 
puted station, x = 0, the v profile was still slightly affected by the inaccurate initial v 
profile. Another parabolic calculation was made with Ax = 0.025. Thus, the x = 0 sta- 
tion was the third station downstream and the v profile was much more accurate than in 
the previous solution. These new u and v profiles at x = 0 shown in figure 3 were 
used as inflow profiles in the calculations described in the present paper. Figure 3 also 
shows profiles from two stations farther downstream which were used as outflow bound- 
aries in the Navier -Stokes calculations. This figure shows that the shear layer spreads 
significantly in the computational domain considered. 

Table I summarizes the results of calculations from case 1 made with the new para- 
bolic code profiles. The grid was 10 points (x-direction) by 122 points (y -direction) 
with uniform spacing Ax = Ay = 0.025; thus, the computational domain extended to 
Xmax ~ 0.225 and to yj^^^x = 3.0125. The computed steady-state u and v profiles at 
X = 0.200 are compared in figure 4(a) with the corresponding calculation made with the 
parabolic code. The circles indicate points from the parabolic calculation for which 
Ay = 0.0125; however, not all the points in the calculation are shown. The u profile is 
virtually identical to the parabolic solution. For this x-station, which is one station from 
the downstream boundary, the maximum difference in v in the mixing region is approx- 
imately 5 percent. The pressure profile for the station shown in figure 4(b) is smooth and 
indicates only a slight deviation (about 1 percent) in the mixing region from the constant 
value (indicated by the dashed line) assumed in the parabolic calculation. Thus, the use of 
the new, more accurate inflow profiles removed the pointwise pressure oscillations present 
in the calculations of reference 1. 
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Effect of Courant num ber.- The time step used in the previously described calcula- 
tion was 0.7 times the maximum step allowed for stability in the explicit hopscotch method, 
that is, Nqq = 0.7 where 


N 


Co 


At 

Ax 

lu| + Ivl + \l2c 


(26) 


Since u varies throughout the domain, the time-step size is actually governed by the 
maximum value of u which is U 2 * For = 1.0 (Courant-Freidrichs-Lewy condi- 
tion) for case 1, At was 0.021. For = 0.7, steady-state convergence was obtained 

in 329 steps. As shown in table I, as was increased to 5.0, the number of steps 

required for convergence dropped to 62. In all these calculations, initial condition IC 1 
was used and At was held constant for all time-marching steps. For = 10.0, how- 

ever, the solution diverged. A possible cause of this divergence was roundoff error from 
the tridiagonal matrix inversion occurring when the coefficient matrix did not possess 
diagonal dominance. The diagonal dominance condition, that is. 



(27) 


is a sufficient but not necessary condition for convergence of the matrix reduction. Keller 
(ref. 16) presents an inductive proof that diagonal dominance insures no error growth in 
the Thomas algorithm. Hirsh and Rudy discuss diagonal dominance further in refer- 
ence 17. In the present = 10.0 calculation, the continuity equation was not diago- 

nally dominant for any row in the horizontal sweep. In addition, the two momentum 
equations lacked diagonal dominance for many rows. 

Even very large time steps were not possible when the initial flow field was taken 
to be the converged steady-state flow field. A calculation was made in which was 

increased by a factor of 1.05 every time step starting with = 1.0. When 

reached 54.6, the solution began to diverge rapidly in the next few time steps. 

Effect of downstre am boundary condition.- To test the effect of the linear extrapola- 
tion outflow boundary conditions, a calculation was made for case 1 with = 10^, 

O' = 0, and = 1.0 with 20 x-stations instead of 10 so that the domain extended to 

X = 0.475. The u and v profiles at x = 0.200 (a station which was now near the 
middle of the domain instead of one column from the boundary) were virtually identical to 
those of the previous calculation. The maximum difference in u was less than 0.5 per- 
cent and for v less than 1.0 percent. Thus, for this problem, linear extrapolation is a 
reasonable boundary condition. 

Comparison with e xplicit method.- The ADI results for case 1 with Nj^g = 10^ 
were compared with solutions of the same equations in the same computational domain 
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obtained by using the explicit hopscotch method described in reference 1. For = 0.7, 

the hopscotch solution required 238 steps to reach steady-state convergence. The com- 
puted steady-state flow fields were virtually identical, with maximum differences of 
much less than 1 percent. For each time step, the CPU time on the Control Data 6600 
computer system for hopscotch was 1.08 x 10"^ sec/node, and for the ADI method, 

3.74 X 10“3 sec/node. (Neither code was optimized.) The best ADI solution ^62 steps for 
Nqq = 5.0^ took about one-fourth as many steps as the hopscotch method, so that the total 
CPU time for the two methods was comparable. (The hopscotch method was unstable for 
Nqq = 0.9 for this problem.j 

Effe ct of initial conditions .- Calculations were also made with initial conditions 
which were not considered to be as good an approximation to the steady-state result as 
was IC 1. The steady -state downstream values of u and v were multiplied by 0.8 in 
IC 2 and 0.6 in IC 3. Linear Interpolation was again used to compute interior values. As 
shown in table I, convergence was actually reached 10 steps sooner with IC 2 and 30 steps 
later with IC 3 than with IC 1. However, in both cases the steady -state flow fields were 
identical to the one computed for IC 1. 

Effect of increased Re ynolds number . - The results of calculations for case 1 with 
NRe = 5.0 X 10^ from reference 1 are shown in figure 5. (It should be noted that in fig. 7 
of ref. 1 these plots were incorrectly captioned as being for subsonic-supersonic mixing.) 
As shown in figure 5(a), with Ax = Ay = 0.025 (10 x 122 grid), the u profile at 

X = 0.15 was accurately predicted, whereas the v profile exhibited an oscillation 
near its maximum value in the viscous mixing region. Halving the grid so that 
Ax = Ay = 0.0125 (19 X 122 grid) eliminated this oscillation. The pressure profiles 

shown in figure 5(b) have small pointwise oscillations for both grids although the maximum 
deviation from the parabolic code value is less for the finer grid. New initial profiles 
were also generated for case 1 (Nr^ = 5.0 x lo3^ by recomputing the parabolic solution 
with a smaller Ax step. As before, the parabolic solution initial profile for u was an 
error -function curve fit between Uj and U2; a zero-v initial profile was used. The 
results from the ADI calculations with these new inputs are summarized in table II. A 
19 X 122 grid was again used with x extending to 0.225 and y extending to 1.50625. A 
comparison of velocity profiles with the parabolic calculation for the mixing region is 
shown in figure 6(a) for x = 0.15. Excellent agreement is again obtained for both u 
and V. As shown in figure 6(b), the new pressure profile is smoother than the reference 1 
profile although very slight oscillations are still present. (The y values of the ref. 1 
profile have been shifted in fig. 6(b) to correspond to those of the new calculation.) The 
deviation from the free-stream value is still much less than 1 percent. 

For N^q = 0.9, the ADI solution required approximately five times as many steps 
to reach steady state as the solution with the hopscotch method. For Nq^ = 0.9, the new 
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inflow profiles did not affect the convergence rate; however, for = 5.0, ADI required 

103 steps (3 times as many as hopscotch) with the new inputs instead of the previous 
134 steps. The solution diverged for = 10.0. This divergence was again possibly 

due to the loss of diagonal dominance in one or more of the coefficient matrices. 

Some tests were also made for case 1 with = 5.0 x 10^ to determine the effect 

of adding artificial viscosity to the continuity equation. In reference 1, it was shown that 
small values of o in equation (17) smoothed the pressure in the hopscotch calculations 
of case 2 at a higher Reynolds number. With a = 0.25 and 0.50, the u and v profiles 
were not altered for = 0.9 while the slight pressure oscillations were smoothed. 

For a = 0.50, the oscillation was reduced to one point. However, as indicated in table n, 
the convergence rate decreased about 10 percent as <y was increased from 0 to 0.5. In 
addition, for = 5.0, a computation with a = 0.5 also diverged as had the computa- 

tion with 0 = 0, For all computations in which artificial viscosity was used, a constant 
value of a was used during the entire calculation. 

Calculations were also made for case 1 with = 8. 1 x 10^, A 19 x 122 grid 

(xmax = 0,225 and ymax = 1.50625) was again used and inflow profiles were generated 
by using the parabolic code. The results are listed in table II and the computed steady- 
state profiles are shown in figure 7 for x = 0.1875, a station near the outflow boundary. 
The agreement with the parabolic solution is very good; for example, the maximum devia- 
tion from the free-stream pressure is less than 0.1 percent. Even though the shear layer 
spreads only slightly over the solution domain considered, this calculation demonstrates 
that the ADI method is stable for significantly high Reynolds number flows. Initial condi- 
tion IC 1 for this flow gives a flow field which closely approximates the steady-state flow 
field; however, over 600 steps were needed to reach steady state with o = 0.25. A solu- 
tion with 0 = 0 was run but convergence at all grid points was not obtained after 
1800 time steps. The calculation was not continued further. The solution diverged for 

Nco = 5-0- 

Effect of diagonal dominance. - Two other forms of spatial differencing which give 
diagonally dominant coefficient matrices were also investigated. These forms were upwind 
differencing and the differencing scheme of Khosla and Rubin. 

Upwind differencing: In upwind differencing, one-sided, rather than centered, finite- 
difference approximations were used for the convective terms in the continuity equation 
and the two momentum equations. For example, the term pvu^ was differenced as 


Thus, the one-sided difference is always taken on the upstream, or upwind, side of the 
point at which the derivative is being computed. This technique has been used extensively 
in previous numerical calculations. Roache (ref. 8) discusses many of these applications 
in detail. The elimination of central differencing for first derivatives gives diagonally 
dominant coefficient matrices in all three difference equations for both sweeps. 

3 

Upwind differencing was used to compute case 1 with = 10 . These results 

are listed in table I. For = 1.0, convergence was obtained in 132 steps, which is 

about one-third as many steps as with central differencing. For = 5.0, only 30 steps 

were needed; this is less than one -half the number required for central differencing 

= 5.0) and one-eighth the number for hopscotch = 0.7). In both computations 

with upwind differencing, the u profiles were virtually identical to those in the central 
differencing calculation, and the v profiles differed by less than 2 percent in the mixing 
region. As with central differencing, the pressure deviates from the free-stream value 
by less than 1 percent. Apparently Ax and Ay were sufficiently small that the first- 
order differencing did not significantly alter the results. Although the numerical diffusion 
introduced by upwind differencing did not affect the profiles, it may have been at least 
partially responsible for the improved convergence rate. 

A converged solution could not be obtained for = 10.0, as was the case with 

central differencing. Oscillations which appeared early in the solution in all variables 
eventually resulted in divergence. The exact cause of these oscillations is not known, but 
it obviously is not the loss of diagonal dominance. 

The results for case 1 with = 5.0 x 10^ are listed in table n. For = 5.0, 

only 48 steps were needed to reach steady state. As with = 10^, the use of upwind 

differencing reduced the total number of steps required by a factor of 2. The u, v, and 
p profiles are again virtually identical to the central differencing results. The artificial 
viscosity inherent in the upwind differencing did not smooth the small pointwise pressure 
oscillations shown in figure 6(b). As with central differencing, the solution diverged for 
Nco = lO.O. 

4 

For the highest Reynolds number, = 8.1 x 10 , the upwind solution converged 

almost four times as fast as the central differencing solution for = 1.0. As shown 

in table n, adding artificial viscosity improved the convergence rate slightly. 

Khosla-Rubin differencing: Another indication that divergence is not always 
caused by a loss of diagonal dominance was given by the results of tests using the differ- 
encing scheme of Khosla and Rubin (ref. 18). This method, which is second-order accu- 
rate and unconditionally diagonally dominant, uses a modified differencing for convective 
terms instead of central differencing. For example, for the term PvUy, the differencing 
becomes 
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This expression reduces to the usual central differencing at steady state when u?t^ = u? .. 
In the present application it was used only for the implicit convective terms in each equa- 
tion; the central differencing form was retained in the Dj and terms. The method 
was shown to be unconditionally stable for a one -dimensional model equation (Burgers’ 
equation) in reference 18 using a linear stability analysis, 

3 

As shown in table I, for case 1 with = 10 , steady state was reached in 

129 steps for = 1; this is approximately the same number of steps required with 

upwind differencing. The u and v profiles were again essentially identical to the cen- 
tral differencing results, and the pressure again deviated from the free-stream value by 
less than 1 percent. For = 5,0, however, oscillations developed in the solution which 

produced a rapid divergence even though diagonal dominance was maintained. The same 
result occurred with = 5,0 x 10^ and = 5,0, Similar behavior was observed 

in solutions of the nonconservative form of the nonlinear Burgers’ equation (ref, 19) but 
not for the linear Burgers' equation even when "wiggles” were present in the solution. 

Thus, these examples suggest that the method may not always be stable for nonlinear 
equations. 


Subsonic -Supersonic Parallel Mixing (Case 2) 

Calculations for case 2 were made with central and upwind differencing with 

Npjg = 10^, The boundary conditions shown in figure 2 were used with necessary specified 

values obtained from the parabolic solution. The grid was 10 x 122 with Ax = Ay = 0.05 

so that x^-^ = 0.450 and = 6.025. Initial condition IC 1 was used for all calcu- 

iiicix nicix 

lations. The results of the calculations are summarized in table in. 

The computed steady-state velocity profiles at x = 0,15 are compared in figure 8(a) 
with the corresponding calculation obtained by using the parabolic code. The u profile is 
virtually identical to the parabolic solution, and the v profile shows very good agreement. 
This agreement is also very good at all other x-stations in the solution domain. The 
static-pressure profiles at two stations, x = 0,15 and 0.30, are shown in figure 8(b). The 
maximum deviation from free-stream static pressure occurs near the sonic point close to 
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the inflow boundary. As shown, the maximum deviation decreases with downstream dis- 
tance. These pressure profiles were smooth with no pointwise oscillations. Gottlieb and 
Gustafsson (ref. 15) observed pressure oscillation in their hopscotch calculation of a shear 
layer flow with a larger transonic region. They were able to smooth the pressure by 
altering the subsonic inflow boundary conditions through specification of appropriate char- 
acteristic variables at the boundaries. Although the pressure does not oscillate spatially 
in the present calculation, incorporation of the Gottlieb -Gustafsson boundary condition may 
improve the pressure near the sonic point in the ADI calculations. 

With central differencing it was necessary to use some artificial viscosity in the 
continuity equation to obtain convergence. With a = 0, the solution for = 0.5 

diverged after 600 time steps; however, with ot = 0.75, the solution converged in 
807 steps. Hopscotch required 683 steps to obtain convergence for = 0.5 and 

cx = 0,25, (Artificial viscosity was also needed with hopscotch.) As shown in table III, 
the number of time steps required decreased as was increased. For = 5,0, 

only 100 time steps were required for convergence. The solution diverged when 

= 10.0. It can also be noted from these results that the presence of subsonic flow 
increased the total number of time steps required in comparison with the fully supersonic 
case 1 results with = 10^. 

ADI with upwind differencing required no artificial viscosity (i.e., a = 0) for 
= 2.5. Steady state was reached in 95 steps, which is about one-seventh of the steps 
required for hopscotch for = 0,5, No solution could be obtained for N^^^ = 4.0 

with O' = 0, 0.25, or 0.75. 


CONCLUDING REMARKS 

A sequential (uncoupled) time-asymptotic finite-difference alternating-direction 
implicit (ADI) method was tested on laminar parallel free mixing flows by using the com- 
pressible Navier-Stokes equations. For the mixing of two supersonic streams with a 
Reynolds number of 10^, a converged steady -state flow field was computed with a Courant 
number of 5.0 and central differencing. The solution diverged for a Courant number of 
10.0, caused perhaps by the loss of diagonal dominance in the coefficient matrices of the 
governing equations. There are other factors, however, which may also have affected 
the size of the maximum allowable time step. For example, the effect of lagging the pres- 
sure and cross -derivative terms is not known. In addition, the nonlinear coefficients in 
the convective terms were lagged in order to linearize the equations. Such a linearization 
may only be appropriate for small time steps. Therefore, a fully coupled procedure with 
a better linearization technique may be needed to gain somewhat larger time steps. 
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Two procedures which give unconditional diagonal dominance were also tested and 
both were found to be limited in allowable time step size. Solutions with upwind differ- 
encing of the convective terms diverged for a Courant number of 10.0; whereas solutions 
with the Khosla-Rubin method diverged for a Courant number of 5.0. In both of these solu- 
tions the failure cannot be attributed to the loss of diagonal dominance. It appears that the 
methods are not unconditionally stable for the full nonlinear Navier-Stokes equations, as 
the linear stability analysis for model equations indicates, but further research is required 
before a definite conclusion can be reached. Perhaps some alternate form of linearization 
may give a more stable scheme even with an uncoupled solution procedure. 

The use of upwind differencing in the ADI method significantly improved the conver- 
gence rate. With a Courant number of 5.0, the upwind version of ADI converged in about 
one-half the number of time steps required with central differencing and one-eighth the 
number required by the explicit hopscotch method. Thus, the expected advantage of an 
implicit method was realized for the supersonic mixing case with a Reynolds number 

O O 

of 10 . When the Reynolds number was raised to 5.0 x 10 , the hopscotch method proved 
to be superior to the ADI in terms of convergence rate and total computing time. 

The ADI method also converged in less total computing time than hopscotch for the 
parallel mixing of a subsonic stream with a supersonic stream for a Reynolds number 
of 10^. Thus, it appears that ADI can be competitive with explicit methods for some prob- 
lems with regard to total computing time; however, the ADI method requires more than 
twice the computer storage that hopscotch requires, and ADI is more complex to code 
than most explicit methods. This result cannot, of course, be generalized until compari- 
sons are made over a wider range of flow problems. 

With central differencing for the convective terms, it was necessary to add artificial 
viscosity to the continuity equation to obtain convergence for subsonic -supersonic mixing. 
None was needed, however, for upwind differencing, which inherently possesses some 
artificial viscosity. For completely supersonic mixing, artificial viscosity was necessary 
only for the highest Reynolds number, 8. 1 x 10'^, when central differencing was used. 

The very large time steps which were expected prior to the present investigation 
were not possible for the uncoupled procedure for the test problems considered. A fully 
coupled solution technique may not yield a significantly larger time step than an uncoupled 
procedure. For nonlinear equations, the high-frequency waves in the solution domain are 
coupled with the low-frequency waves. These high-frequency waves can grow undamped 
and lead to divergence if the solution and/or method parameters are not highly smooth. 

It is possible that too large a time step may adversely alter the wavelength of information 
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coming into the solution domain and cause such a divergence. K this happens, coupling 
the equations will probably not prevent the divergence and allow larger time steps. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
August 3, 1976 
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TABLE L- SUMMARY OF RESULTS FROM CALCULATIONS 


FOR CASE 1 AT = 10^ 


[Ax = Ay = 0.025; = 0.225; grid dimensions, 10 x 12^ 





Time steps required for converge 

^Co 

IC 

a 

ADI 

ADI 

ADI 




central 

upwind 

Khosla-Rubin 

0.7 

1 

0 

329 



1.0 

1 

0 

202 

132 

129 

5.0 

1 

0 

62 

30 

Diverged 

10.0 

1 

0 

Diverged 

Diverged 


.7 

2 

0 

319 



.7 

3 

0 

352 




Hopscotch 
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TABLE II.- SUMMARY OF RESULTS FROM CALCULATIONS 
FOR CASE 1 AT HIGH REYNOLDS NUMBERS 


[ax = Ay = 0.0125; = 0.225; grid dimensions, 19 x 12^ 


^Re 

^Co 



Time steps required for convergence 

IC 

a 

ADI 

central 

ADI 

upwind 

ADI 

Khosla-Rubin 

Hopscotch 

5.0 X lo3 

0.9 

4 

0 

168 



33 


.9 

4 

.25 

175 





.9 

4 

.5 

184 





5.0 

4 

0 

103 

48 

Diverged 



10.0 

4 

0 

Diverged 

Diverged 




10.0 

4 

.5 

; Diverged 




8.1 X 10^ 

1.0 

1 

.25 

666 

183 




5.0 

1 

.25 

Diverged 

Diverged 




. 1.0 

1 

0 


199 
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TABLE m.- SXJMMARY OF RESULTS FROM CALCULATIONS 
FOR CASE 2 AT = 10^ 

[ax = Ay = 0.05; = 0.450; grid dimensions, 10 x 12^ 


Nco 

0.5 

.5 

.5 

.9 

2.5 

4.0 

5.0 

10.0 
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Computational domain 



Case 1 Case 2 

Mj = 1.68 Mj = 0.11 

M 2 = 3.00 Mg = 3.00 


Figure 1.- Standard test problems for mixing of two parallel streams. 



Subsonic inflow 


y = y 


max i 


Subsonic inflow 

u,v specified 
P^ = o 


Supersonic inflow 
p,u,v specified 


y = 0 



X a 0 


X = X 


max 


Line of symmetry 
V = 0 

u = p =0 

y ^y 


Figure 2.- Schematic of computational domain with boundary conditions for case 2 
and grid alinement at upper and lower boundaries. 
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= 0 



u V 

3 

Figure 3.- Velocity profiles for case 1 for = 10 computed with parabolic 

shear layer code of reference 14, 




3.7 3.9 4.1 X 10" 

P 


ADI 

Parabolic code 


(b) Steady-state pressure profile. 
Figure 4.- Concluded. 




ADI 



Ax = Ay = 0.0125 


u V 

(a) Steady-state velocity profiles. 

Q 

Figure 5.- Supersonic-supersonic (case 1) mixing for = 5.0 x 10 . 

X = 0.15. (Calculation from ref. 1.) 
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(b) Steady -state pressure profiles. 
Figure 5.- Concluded. 
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